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ABSTRACT 

The observed giant outbursts of Luminous Blue Variables (LBVs) may 
occur when these massive stars approach their Eddington limits. When this 
happens, they must reach a point where the centrifugal force and the radiative 
acceleration cancel out gravity at the equator. We call this the f2- limit. When 
stars are close to the tt- limit, strong non-spherical mass loss should occur. This 
suggests a scenario where a slow and very dense wind, strongly confined to the 
equatorial plane, is followed by a fast and almost spherical wind. We compute 
two-dimensional hydrodynamic models of the evolution of the nebula formed 
from such interacting winds, using parameters consistent with the outburst 
of rj Carina in the last century. This outburst gave birth to the Homunculus, 
the hourglass-shaped inner part of a highly structured circumstellar nebula. 
Assuming the star was very close to the Q-limit during outburst, our models 
produce gas distributions that strongly resemble the Homunculus on large and 
small scale. This supports the general conjecture that giant outbursts in LBVs 
occur when they approach the Eddington limit. Our models constrains the 
average mass loss rate since the outburst to values smaller than the present-day 
mass loss rate and suggest that 77 Car is approaching another outburst. Our 
models imply that the occurrence of giant LBV outbursts depends on the initial 
stellar rotation rate, and that the initial angular momentum is as important to 
the evolution of very massive stars as their initial mass or metallicity. 

Subject headings: Stars: Circumstellar Matter — Stars: Mass Loss - 
Stars: Rotation — Stars: Evolution — Stars: Individual: rj Carinae 
— Hydrodynamics 
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1. Introduction 

Current stellar evolution models predict that Galactic stars initially more massive than 
25-30M© lose more than 50% of their initial mass, and stars above 3O-35M more than 
80% (Maeder 1992, Woosley, Langer & Weaver 1993). Much of this mass loss is thought 
to occur during a short-lived and highly unstable stage preceding the Wolf-Rayet stage, 
that may be observed in the form of the luminous blue variables (LBVs; cf. Maeder 1989, 
Pasquali et al. 1997), which are located close to the Humphreys-Davidson (HD) limit in the 
HR diagram, beyond which no normal stars are observed (Humphreys & Davidson 1979). It 
is strongly debated what produces the giant LBV outbursts observed in these stars (Langer 
et al. 1994, Nota & Lamers 1997). Among the potential mechanisms (see Humphreys & 
Davidson 1994, and references therein), the idea that massive stars reach their Eddington 
luminosity close to the HD limit is particularly appealing since the Eddington limit appears 
to be located very close to the HD limit in the HR diagram (cf. Davidson 1971, Lamers & 
Noordhoek 1992). 

Time-varying massive star winds are able to produce circumstellar nebulae with 
kinematic properties similar to those observed around LBVs (Frank, Balick & Davidson 
1995; Nota et al. 1995; Garcia-Segura, Mac Low & Langer 1996, Frank, Ryu, & Davidson 
1998). The morphology of the nebula around the extraordinary star ^Carinae may 
give essential clues for the understanding of these giant outbursts. During its eruption 
from 1840 to 1860 A.D., r^Car — today a telescopic object — was the second brightest 
star in the sky (van Genderen & The 1984). Recent observations by the Hubble Space 
Telescope have revealed in spectacular detail the resulting circumstellar nebula (Humphreys 
& Davidson 1994, Morse et al. 1998), the hourglass-shaped inner part of which is known as 
the Homunculus (Meaburn, Wolstencroft & Walsh 1987, Allen 1989). This is now the best 
studied example of the bipolar structures often observed around LBVs such as rj Car (Nota 
et al. 1995). 
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2. Winds of rotating stars near the Eddington limit 

It has been proposed (e.g. Maeder 1989) that sufficiently luminous stars, after core 
hydrogen exhaustion, may arrive at or exceed their Eddington limit 

r = L/L Edd = i, (i) 

as their Eddington luminosity L Edd = 4ncGM/ k drops below their actual luminosity L as 
the opacity coefficient k increases. The theory of radiation-driven stellar winds predicts 
that the mass loss rate will increase as M oc (1 — with fi > for stars approaching 
the Eddington limit (cf. Castor, Abbott & Klein 1975). This mass loss rate formally 
diverges for T — > 1 (though see Owocki & Gayley 1997), suggesting that the strong mass 
loss associated with T ~ 1 may be related to giant outbursts of LB Vs. 

However, rotation reduces the luminosity required for all external forces to balance 
each other at the stellar surface (Langer 1997). The Eddington limit T < 1 should then be 
replaced by a criterion that we will call the fi-limit, 

^ = v rot /v crit < 1, (2) 

with i£ it = Ve SC /2 = GM(1 — T)/R, M and R being the stellar mass and radius, and 
f esc being the polar escape velocity. Q = 1 implies that centrifugal and radiation force 
balance gravity at the equator, while at higher latitudes gravity still dominates. The 
feedback of rotation on the local surface luminosity is neglected here, since according to 
the generalized von Zeipel theorem (Kippenhahn 1977) the radiation flux at the equator 
may be either reduced or enhanced, depending on the internal rotation law; this may have 
been overlooked in the recent criticism of the fi-limit by Glatzel (1998). Here, we apply the 
result of Friend & Abbott (1986) that the mass- loss rate of rotating hot stars depends on f2 
as M oc (1 - n)~ v with v ~ 0.43. 

We divide the evolution of a star approaching the f2-limit — that is, going through an 
outburst — into three phases. In the first phase, before the star reaches the fi-limit (Q < 1), 
it has the fast, energetic wind expected of a luminous blue star. In the second phase, 
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reaching the fi-limit (Q ~ 1) has three consequences for the wind: the mass loss rate is 
much higher than before; the mass flux increases strongly at latitudes close to the equator; 
and the bulk of the wind is slow since the equatorial escape velocity is almost zero. In the 
third phase the outburst is over, the stellar radius has decreased, and the configuration is 
similar to the first phase. The smaller value of T before and after the outburst has two 
consequences: larger wind velocities due to larger escape velocities and smaller values of fl 
leading to more spherical winds. 

To compute the latitudinal dependence of the wind properties of a star close to 
critical rotation ideally requires mult i- dimensional models of the star and its outflowing 
atmosphere, which are not available. However, Langer (1997, 1998) argued that the stellar 
flux and the radius might still vary only weakly from pole to the equator in very luminous 
stars. Therefore, we applied equations similar to those found by Bjorkman & Cassinelli 
(1993, BC) for winds of rotating stars in the limit of large distance from the star: 

v 00 (6) = (v csc (i-nsmey , (3) 

(47rr 2 p)oc(#) = !<JM (1 -SI sm9)* / Voo (9) , (4) 

where we set the parameters defined in BC to ( — 1 ,7 = 0.35, and £ = —0.43. The 
correction factor 5 is introduced to ensure that the total stellar mass loss rate M obeys 
Mo = J v OQ (9)p(9) sin9 d9d(p at the inner boundary of our grid (cf. Table 1). is the 
terminal wind velocity, and (47rr 2 p) QO the terminal wind density times Anr 2 , as function of 
the polar angle 9. The quantity a is defined by 

a = (cos (j)' + cot 2 9 (l + 7 - ^ sinfl \ , ^ A ^ 
\ \ 1 — il smV J J 

with (/>' — Q sm9v crit /(2\ / 2v 00 (9)). Eq. (5) differs from the corresponding quantity defined 
by BC in their implicit formula (26) (cf. also Owocki et al. 1994, Ignace et al. 1996). 
This difference came along originally through a misinterpretation of BC's equations, i.e., in 
equations (3) to (5) and in the equation defining 0', the #s were taken as 9 s, the initial 
co-latitude of the streamline. We note that since for wind compressed zone models near 
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the equator it is 6 ~ 9, our models are similar to wind compressed zone models when 
Q ~ Vl th ~ 1, where f2 t h is the threshold value for the formation of a wind-compressed disk. 
Since <fi' < n/2 for Q < 0.995, our formulation has the advantage of avoiding the formation 
of wind compressed disks for large Q, a structure which can not be numerically resolved 
in our calculations, whose properties cannot be well predicted, and whose very existence 
has even been questioned (Owocki et al.1994, 1996). At the same time, wind density and 
velocity distributions obtained from our approach are similar to those derived from the 
formalism of BC, provided that fl t h ~ 1 and we choose Qrc = ^^th, where Qbc is the 
value to be inserted in BC's equations. We shall see that the exact nature of the latitude 
dependence of the wind properties is not essential for our main results, as long as a dense 
wind with enhanced mass loss rate close to the equator occurs between two phases of an 
energetic, more or less spherical wind. 

We simulate the LBV outburst phenomenon by assuming the wind properties to be 
constant during each phase. For the pre-outburst wind, which is only used to initialize 
the numerical grid and to which our results are insensitive, we took a mass-loss rate 
of M = lO _3 M yr _1 , a wind final velocity = 450 kins" 1 and Q = 0.53. For 
the post-outburst wind, we used two different sets of parameters. The one that we 
prefer appears to reasonably reproduce the observed morphology of the Homunculus, 
with M — 1.7 x lO~ 4 M yr -1 , = 1800 kms" 1 , and Q = 0.13. For comparison we 
also computed a model that represents rj Car's presently observed wind parameters of 
M = 3 x lO^Moyr" 1 , Voo = 800 km s _1 , and Q = 0.3 (Davidson et al. 1995). The high 
wind velocity of our preferred model, = 1800 kms -1 , corresponds to a stellar radius 
R ~ 21R Q or logT cff ~ 4.7 for an O star wind with ( = 3 at log L/L Q = 6.4, implying that 
the star strongly contracted after the episode of mass-loss. 

The parameters of the outburst wind largely determine the morphology of the resulting 
nebula. To compute these parameters, we assumed the following stellar properties: 
M = 8OM , log L/Lq = 6.4, and logT c fj = 4.2, implying R = 21Oi? , in agreement with 
observational estimates (cf. Humphreys & Davidson 1994). The final wind velocity follows 
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from the escape velocity for a specified value of T. For our preferred model we took a mass 
loss rate of M = 7 x lO _3 M yr _1 and obtained a Homunculus mass of roughly O.15M 
(van Genderen & The 1984), while for our model assuming present-day wind parameters, 
we used an increased outburst mass loss rate of 5 x 1O~ 2 M yr" 1 such that we obtained a 
nebula mass of 1 M (Humphreys & Davidson 1994). 

3. Hydrodynamic models 

We perform two-dimensional hydrodynamic simulations of the wind interaction; first 
results have already been reported by Garcia-Segura, Langer & Mac Low (1997). We use the 
hydrocode ZEUS-3D developed by M. L. Norman and the Laboratory for Computational 
Astrophysics. ZEUS-3D is a finite-difference, fully explicit, Eulerian code descended from 
the code described by Stone & Norman (1992). We used spherical coordinates for our 
simulations, with a symmetry axis at the pole, and reflecting boundary conditions at the 
equator and the polar axis. See Garcia-Segura et al. (1996) for further details about our 
numerical method. Our models have grids of 200 x 360 zones, with a radial extent of 0.125 
pc, and an angular extent of 90°. The innermost radial zone lies at r = 9.7 x 10 15 cm. 

We compute the hydrodynamic evolution of the circumstellar gas, starting our 
computations at a time t = 1840 yr, and run them until t = 1995 yr. Our outburst scenario 
leads to a characteristic distribution of the circumstellar gas. The initial fast wind blows 
a stellar wind bubble which forms the background for the subsequent development of 
the nebula. During outburst, the wind becomes slow and dense, and the stellar rotation 
concentrates it toward the equatorial plane (BC; Ignace et al. 1996). When the final fast 
wind starts in the center of this nebula, it sweeps up the dense wind from the outburst into 
a thin, radiatively cooled shell that fragments due to dynamical instabilities (Garcia-Segura 
et al. 1996). The shell expands more easily into the lower density wind at the poles, 
producing a double-lobed structure, as shown in Figures 1 and 2. 

In Figure 1, we show six models computed with various values of Q and T for the 
outburst wind given in Table 1, and otherwise using the parameters of our preferred model. 
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We find three major results: First, the nebular shape appears nearly independent of the 
Eddington factor, T. Second, the nebula is strongly confined in the equatorial plane only 
for the case with nearly critical rotation (fl = 0.98). Finally, in this case, we obtain a 
structure very similar to that of the Homunculus, with two almost spherical lobes and an 
equatorial density enhancement that is an expanded relic of the outburst wind. In Figure 2 
we show the results of a computation identical to that of Figure la — in particular with 
Q = 0.98 — but using twice the radial resolution. We find it striking how well this model, 
without much fine tuning, not only reproduces the large-scale, bipolar morphology, but also 
the small-scale turbulent structure seen in high-resolution observations of the Homunculus 
(Humphreys & Davidson 1994, Morse et al. 1998). 

We also compute a model at the same resolution as Figure 2 using r\ Car's presently 
observed wind parameters. We find that the large scale shape of the resulting nebula, 
shown in Figure 3, is almost identical to that shown in Figure 2, but the higher wind 
densities cause the wind termination shock to be strongly radiative. This changes the 
Vishniac instabilities (Vishniac 1983) seen in Figure 2 into ram-ram-pressure instabilities 
(Vishniac 1994, Garcia-Segura et al. 1996), which have a much spikier morphology and a 
shorter wavelength, inconsistent with the observed structures. We emphasize that only 
the properties of the post-eruption wind are responsible for this feature, not the larger 
shell mass obtained in this case. This result might imply that the current wind is not 
representative of the wind over the last 140 years. Instead, 77 Car might have been smaller 
and hotter in the recent past, with a post-outburst wind that has become slower and more 
dense with time. This is consistent with its gradual visual brightening over the last 140 
years (Humphreys & Davidson 1994), and suggests that it is evolving towards another giant 
eruption. 

4. Discussion 

Our work extends previous hydrodynamic models for bipolar LBV nebulae (Nota et al. 
1995, Frank et al. 1995, 1998, Dwarkadas & Balick 1998) by relating the outburst to the 
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properties of evolving, massive, post-main sequence stars. In contrast to Nota et al. (1995) 
and Frank et al. (1995), who concluded that a strong equatorial density enhancement 
must have existed before the outburst occurred, we obtain the two lobes, including 
their small-scale structure, and the equatorial density enhancement self-consistently as a 
consequence of the evolutionary state of the star. Frank et al. (1998) used an arbitrary 
non-spherical wind during the post-outburst phase to produce a bipolar nebula. However, 
such a wind will only produce a bipolar shape if the wind termination shock is strongly 
radiative and therefore momentum conserving, a condition we have shown to be inconsistent 
with the small-scale morphology of the Homunculus. Dwarkadas & Balick (1998) introduce 
instead a ring-like density distribution, again without relating it to the underlying star. 

Our result appears to be quite general, because all stars, even slow rotators, must by 
definition arrive at critical rotation if they approach their Eddington limit. The strong 
dependence of the nebula shape on Q shown in Figure 1, as well as the clear bipolar nature 
of virtually all LBV nebulae (Nota et al. 1995), lends strong support to the idea that LBV's 
are stars approaching their Eddington limits that reach critical rotation and lose large 
amounts of mass quickly. In fact, a general mechanism for giant LBV outbursts is needed 
in order to understand the absence of stars beyond the Humphreys-Davidson limit in the 
HR diagram, and the bipolar nature of most LBV nebulae. Therefore, even though bipolar 
nebulae may also form from interacting binary stars (e.g. Han, Podsiadlowski & Eggleton 
1995), and binarity has been repeatedly proposed also for 7] Car (van Genderen, de Groot 
& The 1994, Damineli, Conti & Lopes 1997), it appears useful to continue to pursue single 
star models. The conjecture that the Homunculus nebula around r\ Car is a paradigm rather 
than a freak has recently been supported by its strong similarity to the nebula around the 
LBV HR Carinae, as found by Weis et al. (1997) and Nota et al. (1997). 

NL is grateful to Gloria Koenigsberger for enabling an extended visit to the UNAM 
Astronomy Institute, and to many colleagues there for their extremely warm hospitality. 
We thank J. Fliegner, S. Owocki, M. Peimbert, S. White and in particular J. Bjorkman for 
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Figure Captions 

Fig. 1. — Logarithm of the circumstellar gas density (in g cm -3 ) in our hydrodynamic models 
at t — 1995 yr. Figs, la-f correspond to models a-f in Table 1. 

Fig. 2. — Same as Figure la, but using twice the radial resolution. The resemblance of 
the Vishniac instabilities in the swept-up lobes with the corresponding features in the high- 
resolution observations by the Hubble Space Telescope (Humphreys & Davidson 1994, Morse 
et al. 1998) is striking. The mass of the Homunculus nebula in this model is ~ O.15M (van 
Genderen & The 1984). 

Fig. 3. — Same as Figure la, but using r\ Car's presently observed wind parameters. The 
higher mass- loss rate changes the nature of the shell instabilities, as discussed in the text, 
producing a morphology that disagrees with the observations. The mass of the Homunculus 
nebula in this model is ~ 1M Q (Humphreys & Davidson 1994). 
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Table 1: Parameters during the outburst phase for the models shown in Figure 1. i> esc , the 
escape velocity from the stellar surface in polar direction, is also the final wind velocity in 
polar direction. See text and figure captions for definitions of other quantities. 
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0.98 


0.50 


186.7 


190.5 


269 
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0.90 


0.50 
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190.5 
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0.98 


0.75 


132.0 


134.7 
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0.90 


0.75 


121.2 


134.7 


190 
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0.98 


0.98 


37.3 


38.1 


54 
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0.90 


0.98 


34.3 


38.1 


54 
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